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Abstract. We study the stochastic dynamics of Ising spin models with random 
bonds, interacting on finitely connected Poissonnian random graphs. We use 
the dynamical replica method to derive closed dynamical equations for the joint 
spin-field probability distribution, and solve these within the replica symmetry 
ansatz. Although the theory is developed in a general setting, with a view to 
future applications in various other fields, in this paper we apply it mainly to 
the dynamics of the Glauber algorithm (extended with cooling schedules) when 
running on the so-called vertex cover optimization problem. Our theoretical 
predictions are tested against both Monte Carlo simulations and known results 
from equilibrium studies. In contrast to previous dynamical analyses based on 
deriving closed equations for only a small numbers of scalar order parameters, the 
agreement between theory and experiment in the present study is nearly perfect. 



1. Introduction 

The interest in studying finitely connected (FC) spin systems on random graphs, as 
introduced in pQ more than twenty years ago, has grown in recent years. For this 
there appear to be at least two reasons. Firstly, FC spin systems can be seen as 
an intermediate step between fully connected mean-field spin models [2] and finite- 
dimensional spin models. Although they are still of the mean field type in the 
sense that random site permutations are irrelevant in the mean field limit, the finite 
connectivity introduces notions of site neighborhood, distances, etc. This attractive 
property has drawn many into this field, and we have by now achieved a thorough 
understanding of the equilibrium behavior of FC spin systems O [4j [5j [6j [7] . Secondly, 
many optimization and decision problems in theoretical computer science can be 
mapped into models of FC spin systems. This mapping allowed such problems to 
be studied with analytic methods of statistical mechanics, and has been very fruitful 
especially in the study of K-SAT [3|9], vertex covering [10], and graph coloring [TT1[T2] . 

Although our understanding of the equilibrium properties of FC spin systems is 
now quite advanced, that of the non-equilibrium behaviour of such systems is, despite 
recent progress [HI HU [T5] [16j H3 OS US], still relatively limited in comparison. 
Dynamical studies are generally harder, by definition, as they incorporate the 
equilibrium state as a special case. In the domain of the dynamics of FC spin systems 
the generating functional method (or path integration technique) of [20] is the only 
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exact method available today. There has been some success in applying this method 
to finitely connected soft spin systems [T31Q3] and Ising spin systems [IB]. However, 
the generating functional method leads one in FC systems to a formalism involving a 
rather complicated dynamical order parameter (describing the joint statistics of single- 
site spin 'paths' and single site field perturbation 'paths') which is generally difficult to 
handle. Even for parallel dynamics [16j it is effectively equivalent to having a number 
of scalar order parameters that grows exponentially with the number of discrete time 
steps considered. For that reason, even in generating functional analysis studies one 
is in practice forced to make further approximations to tame this explosion of order 
parameters. 

An alternative approach to the dynamics of FC spin systems is dynamic replica 
theory (DRT) (2UE2], which was initially developed for fully connected systems. In 
contrast to generating functional analysis, DRT in its present form is not (yet) exact; 
however, one can increase its accuracy systematically by increasing the size of the 
chosen order parameter set [35] . The great advantage of DRT in the study of FC spin 
systems, compared to generating functional analysis, is that the effective number of 
order parameters does not grow with time. Recently, the DRT method [17] and its 
equivalent [15] were used to study the dynamics of FC Ising spin systems, but only 
for a relatively small number of dynamic order parameters. Although its performance 
on regular random graphs was found to be very good [171 115j , for random Poissonian 
graphs it was found to be quite poor [17]. In the present paper we develop the DRT 
method further, and cure the previous limitations by increasing the size of the order 
parameter set, following [22J, [15], to the full joint spin-field distribution. We then 
demonstrate the performance of the resulting improved theory by application to the 
so-called minimal vertex cover problem |10j on Poissonnian random graphs. 

This paper is organized as follows. In section [2] we define our model and derive 
an exact dynamical equation for the joint spin-field probability distribution. In the 
next section [3] we close this equation using the standard assumptions and procedures 
of DRT. We simplify our dynamical theory by making the standard replica symmetry 
ansatz in section T3. 21 In section 3] we apply our resulting formalism to the dynamics 
of the Glauber algorithm, extended with simulated annealing type cooling schedules, 
when running on the minimal vertex cover problem. The outcome of solving our 
dynamical equations numerically arc compared to measurements taken in Monte Carlo 
simulations. Finally, in section [6J we summarize and discuss our results. 

2. Model definitions and macroscopic laws 

We consider a system of N Ising spins, Oj G { — 1,1}, which are placed on the vertices 
of a random Erdos-Renyi graph |23| . Spins interact only when they are connected. 
Their microscopic dynamics are governed by a Glauber type stochastic algorithm. At 
each iteration of this algorithm a site i is drawn randomly from the set {1, . . . , TV} of 
all sites, and spin Oi is subsequently flipped with probability 



o~i ) = ~[1 - o-it&nh[[3hi(a)]] 



(1) 



where hi is a local field, defined as 




(2) 
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with er = (cr 1; . . . , er/v)- The inverse temperature (3 — T _1 controls the level of noise in 
the system; the dynamics is random for (3 = 0, and fully deterministic for (3 — > oo. The 
parameter 9 defines a uniform external field. The set of random variables {cy Jy} is 
regarded as a quenched disorder. The bonds Jy are symmetric, viz. Jy = Jjj, 
and drawn independently from a probability distribution P{J). The independently 
distributed random variables ctj € {0,1} are the entries of a symmetric adjacency 
matrix with zeroes on the main diagonal, defining the random graph. In this paper 
we consider finitely connected (FC) random graphs of the Erdos-Renyi [23] type, where 

Vi<j; P(c ij ) = ^6 Cij , 1 + {l-^)5 Cij , Q (3) 

with c = 0(N°). In the N — >• oo limit the average number of connections per spin 
(or vertex) remains finite, and the distribution of connectivities (or vertex degrees) is 
given by a Poisson distribution with mean c: 

P c (k) = c k e~ c /k\ (4) 

The process (JlJ can be written in the form of a master equation for the evolution of 
the microscopic state probability in continuous time[J 

d N 

foPt^) = ^2\pt(Fitr)wi(Fitr) - pt(tr)wi{(r)] (5) 

i—l 

in which Fi is a spin-flip operator Fitt(cr) = Q(ai, . . . , —<Ji, ■ ■ . , a\/v) and the quantities 
Wi(cr) are the transition rates given by 

Wi(cr) = -[1 - ^tanh^er)]] (6) 

This process evolves towards equilibrium Boltzmann probability distribution Poo(c) ~ 
exp[— (3H(cr)], with the Hamiltonian 

H{(t) = — aiCijJijaj — Oy^o-j (7) 

i<j i 

In general it is not possible to solve the 2 N coupled equations ([5]) directly. Therefore, 
instead of following the evolution of the microscopic distribution pt{cr), one turns to 
alternative descriptions of the dynamics in terms of macroscopic observables. 

For the reasons given in the introduction, we now follow the steps of dynamic 
replica theory [22j . and consider the evolution in time of an arbitrary set of £ 
macroscopic observables Jl(<r) = (^(er), . . . , f^(<x)), where each individual O/^tr) is 
taken to be of order O(N ). We derive a Kramers- Moyal expansion for the associated 
macroscopic probability distribution Pt(fl) = J2cr ~ ^( <T )]Pt( cr ) by inserting the 
master equation §5§ into the time derivative of P t (f2), and expanding the result in 
powers of the 'discrete derivatives' A''(cr) = il^(Ficr) — fl^(cr). This gives: 

fj.=l M i 

2 ^ 



(J,,V—1 

t/3 A3 



+ 0{Nt A d ) (8) 

| This involves formally the introduction of random durations for the individual spin updates which 
are N~ 1 on average, and which for finite N are drawn from a specific distribution 1241 . 
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where we used the sub-shell (or conditional) average 

If the diffusion term in the expansion 1(5)) vanishes for N — > oo, then 1(5)) acquires the 
Liouville form, the solution of which describes the following deterministic flow: 

d n = ( £ ^(<r) [n(F <0 -) - n(<r)] ) n>t (10) 



df 



This is then exact for N — > cxd, but not necessarily closed, due to the presence of the 
microscopic probability pt(tr) in (J9j) . In DRT, in order to close equation (flQ|) . one 
assumes equi-partitioning of probability within the $7 sub-shells, i.e. one takes pt(tr) 
to depend on tj only through fi(<x). The impact of this assumption on the accuracy 
of the theory depends critically on the choice of observables Sl(tr). 

In this paper, our choice of observables f2(er) is, as in [22], the (infinite 
dimensional) set given by the joint spin-field distribution: 

D(s,h;a) = ^Y i S.^S[h-h i {a)] C 11 ) 

i 

We assume that this distribution l|lip is well behaved in the sense that it can be 
evaluated first for a finite number £ of field arguments hp, and that the limit £ — > cxd 
can be taken after the thermodynamic limit N — > oo. For now on we thus have 2£ 
observables D(s, hp; tr) with /j, = 1, . . . ,£ and s £ { — 1, 1}- In order to compute (JTOJ) 
we must work out the discrete derivatives A^(<r) = D(s, hp] FiCr) — D(s, hp] a): 



= ^E^w - h ^ F ^ ^Ev^i^ - (12) 

3 3 

Jj ^2 S sM 3 Cij^S ai:1 S [hp- hj(cr)+ 2Jy] + [fy,- hj(tr)- 2Jij] 

[hp- hj(cr)} | + — {5 s ,_ CTj - <y s>CTi }<5 /li(cr)] 



- 5 [h a - 



Thus A^(er) = 0(N" 1 ), so for iV — > cxd the diffusion term in ([5]) vanishes and 
the macroscopic observables D(s,h u ;ar) evolve deterministically according to ITOl) . 
Inserting (JT2J) into (fTU|) gives us a diffusion equation for the joint spin-field distribution: 

3 1 1 

— D(s, hp) = -[l + s tanh^]] D{-s, hp) --[1-s tanh[/3/i J] D(s, /i M ) 



ij^y dh'[l-s'tanh[/3/i'] 



^■T7 ^ ^',<T i ^,a i C i jg[fe / -fei(<T)](5[/i, A1 -fe i (o-) + 2J tf 



£>:/, 



5^/ d/i'[l-s'tanh[/3/i']] 

^ E WA<^^-^( ff )]%/*-fy( ')]) ( ( 13 ) 



with the sub-shell average 



EqMoQ/W n sM S [Djs, hp) - Djs, h a - a) 
£<r EU ^ [£>( a , ^) - £>(*, hp; a') 
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The non-trivial objects in (fT5|) are the two averages, with angular brackets. To 
compute these efficiently we introduce the following kernel, where s £ {0, s'}, 

A[s,s';h,h';s] = (^^S s , i ^ Staj c ij S[h'-hi(a)]S[h-hj(a)+2Jij^\ (15) 



\cN ' . . - '/ D;t 

y 

For s = the kernel (fT5|) defines the joint spin-field probability of connected sites (a 
similar object was used to study the dynamics of the Ising ferromagnet on a regular 
random graph [15 ). In the limit N — > oo, definition (|15[) allows us to write (fT3"|) as 

5 1 1 

— D(s, fc) = - [1 + s tanh[/?/i]] D(-s, h)--[l-s taahf/3/i]] D(s, /i) 
+ \cYl J dh '^ ~ s'tanh[f3h']]A[s,s';h,h';s'] 

s' 

~\cJ2 J dh'[l- s't&nh[(3ti}]A[s,s'-h,ti;0}. (16) 

s' 

This dynamical equation (fH))) is exact for large N, but not yet closed. Closure requires 
eliminating p t (cr) from (|15p . 



3. Replica analysis of the dynamics 

3.1. Closure and disorder averaging 

To evaluate the right-hand side of (|16p we make the usual assumptions of the dynamic 
replica method. The observables D(s,h^;a) are taken to be self-averaging with respect 
to the disorder at any time, i.e. to depend only on the statistics of the {cijJij} 
rather than their realization. Second, we assume equi-partitioning of the microscopic 
probability within the D(s, h^; cr) sub-shells of the conditional average (fT4|) . These 
assumptions and the equivalence of sites after disorder averaging, lead us to 



A[s, s ; h, h ; s] = lim 



N-l I EcrUr^iD^h^-D^h^cr)} 

N^c C \J2cr'I\r^[D(r,h^~D(T,h^cr')} 

x 6 s ,, ai S s , a2 c 12 S[ti- hi(a)]S[h - h 2 {<r) + 2J 12 S}) (17) 
We eliminate the fraction from the above expression via the replica identity 

Sg^-JSsE-E^n^o as) 

L^CT \ J en (T n a=l 

which leads to 

A\s,s';h,h';s] = lim lim / . . . 

cri cr" 

x 6 s ,^6 s ^c 12 6[h - h^Mh - h 2 {(T l ) + 2J 12 s] (19) 

n 

<X=1 TfJ, I 

We can remove the disorder dependent local fields {hi(a a )} from inside the delta 
functions by inserting into (|19l) the following integral representation of unity: 1 = 
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I Yiai dH°6[H" — hi(<j a )]. Writing the latter delta functions in integral form gives 

N-l 



A[s, s'; h, h'\ si = lim Urn V . . . V [ TT [dfff dhf explihfH. 

N^oo n— +0 C ^ — ' ^ — ' / H L 



CT 1 O" ai 



n n 5 m - ^ e - ^ 



d — 1 



(20) 



x 5 sla] 5 sah 5[h'~Hl\ (c 12 5[h-H 1 2 +2J 12 S] C -^ h " h ^ (Ta A 
After the average over the disorder is taken (see | Appendix A| for details), we then find 



A[s, s'; h, ti; §] = Hrn^ Hm E • • ■ E "Wi^-l / II [ 



2tt 



(21) 



6[ti- H\] Y[s[d(t, M - ^E ^^-•H? 

T jJjOL i 

e-^d'KO J dJ p(j)S[h-H^ + 2JS}e- u ^^^ +il >° ] 
±- £ ( /dJ P(J)e-"S-l« , -3 , +*y-f] - l) + 0(1) 



x exp 



The 0(1) term in the exponent of the last line is independent of {s, s' , h, h' ,§}, and 
can always be recovered from the normalization ^ s , J dhdh'A[s, s'; h, ti] s] = 1. Next 
we achieve factorization over sites in (PHI) upon isolating the density 



P(a, h; {(Ti}, {hi}) = — ^ 5cr,a z S[ h - hi] 

i 

via insertion into (|21[) of the ^-functional unity representation 



(22) 



(23) 



ah 



which gives, with the short-hands (g(J))j — JdJ P(J)g(J) and x • y = J2 a xa y a : > 



A[s, s'; h, ti; s] = lim lim 

N— >oo n— >0 



n[ 



2-k/N 



n[ 

ah 



dP(er, ti)dP(a, h) 



2tt/N 



x e 
x 



exp {iV [i £ D a (r, ^)I»(r, M + i E H"> h ) + °(^) 

T » a ah 

+ \ c E / d/ld/l ' P ( ff > &)*V> h')^ e - iJ ^- ff '+^' <7 l - l) J } 

<T<T' ^ J 

E-E/n[^^- iff -i 

iE X( ,„c (T,^)E, «r,»f*[^-fl' 1 a ]-iE er ^^(o'»^)E < *o-,o- 4 *[ft-ftil 
ai 5. „i8[ti - Hi]/^ - Hi + 2 Js]e- iJ ^ +h ^A (24) 



where a = (<ti, . . . a n ), ai = (a}, . . . a") and similarly for the replicated vectors h, etc. 
We rescale the conjugate integration variables according to P(a, h) — > dh P(a, h) and 
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Z)q,(t, h^) — > Ah^Dair, h^). This converts the sums over h and [i in the exponent of 
(124| into well-defined integrals when dh — ► and I — + oo. We write the resulting path 
integral measure as {dPdPdl?}. Next we define an effective single-site measure M; 

/ h A J,JdHdi 1 Mp,MW/[H,M] 

\ im J2<Tl dHdh M[H,h,<r\0] 

M[H,h,tr\9] = e ^-l H - e ]- i ^a,^D c ,( s ,h^)6 a ^ a sih^-H a ]-iP(cr,h) 

and the function 

*[{P,P,D}] =iJ2^D a (s,h^D(s,h li )+iJ2 /dh P(a,h)P(a,h) 

SpLCX CT 

+ log^ IdHdh M[H,h,tr\0] (26) 
+ \cY, j dhdh ' P (°"' A)(e- iJ ^- <T '+ /l '- <7 ] - l' 



Using these definitions and changing the order of the limits N — > oo and n — > allows 
us to write ([24]) in the form 



A[s,s';M';s] = hm lim [{dPdPdD} e ™[{P,P,f>}]+o(i) (2 7) 



J, MM' 

Finally, with the help of the normalization identity J~)„i Jdhdh'A[s,s';h,h';s] = 1, 
we compute (|27| by steepest descent: 

(5 s>1 5 s , CTl 5[/ l '-F 1 ]5[/ l -^+2Js]e-^^ CT '+^'°-]> JiM ^ 
A|s, s ;n, h ; s\ = lim : — 

/p-iJ[h-<7'+/l-<T]\ 

\ I J, MM , v 

(28) 

where {P, P, I?} are determined by extremization of '5. The functional variation of \& 
with respect to P(<x, h), P(<x, h) and D a (s, h^) leads to the stationarity conditions 

D(s,h) = (6 a ,„J[h-H a ]) M (29) 

P(a,h) = (5 a ,(T'S[h-h'}) M (30) 

P(<x, h) = icJ2 [dh P{ct', ti)(e-u\h-<T'+h'-<T] _ x)j (31) 

CT' J 

The conjugate order parameters D a (s,h) and P(cr,h) are seen to play the role of 
Lagrange multipliers, ensuring normalization of D(s,h) and P(cr,h). The physical 
meaning of the density P(er, h) is not yet clear, due to the presence of the vector h. 

We use equation (|31| to eliminate the conjugate order parameters P(tr, h) from 
the measure M. We assume that D{s, h) is sufficiently smooth in h, such that 
£ M Ah^D a (s, h^)f{h^) -> fdff £» a (s, ff)/(Jf) for * -» °°- This lead s to 

M[H,h,tr\6] = expjih- [H - 0] - D a (a a , H a ) 

a 

+ cJ2 Uh' P{a',h!){e- iJ[h ' T ' +h '' T] -l)j} (32) 

CT' J 
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in the definition of the measure M (|25|) . The replica method requires finally that we 
take the n — > limit in equations (|28M30p . To do this we need to make appropriate 
ansdtze for the density P(<x, /i) and for the conjugate order parameters D a (s,H). 



3.2. Replica symmetry 

We evaluate (|28|) - (j30| upon assuming ergodicity, which translates mathematically into 
the so-called replica-symmetry (RS) ansatz. Firstly, the order parameters D a (s,H) 
depend only on a single replica index and are expected to be imaginary, so we put 

D a (s,H) = ilog d(s,H) (33) 

Second, the density P(<x, h) depends on a discrete and continuous vector in replica 
space. The RS ansatz demands its invariance under any joint permutation of their 
indices, which implies [7] that it must be of the general form 

/n 
{dP} W[{P}]Y[P(<r a ,h a ) (34) 

a=l 

where VT[{P}] is a normalized functional distribution, i.e. J{dP} W[{P}} = 1. The 
RS ansatz (|33I34[) , via its implications for the effective measure (|25[). will enable us 
to take the replica limit n — > in equations ()28H30j) . We insert (|33l34p into (j32|) and 
subsequently expand the exponential function containing Prs(<t, h), leading to 

k f ^ 

M RS [H,h,*\6] = \f~ c I II {dJeP{Ji){dP e }W[{Pi}}} (35) 



fe>0 i=l 
n k „ 

x H {al(a a ,H a )e ih ^ H ^l[[j2 / dhfP t (af, hf)e~ u ^" +h >" 

a=l l=\ a a 



We write averages with respect to the RS measure (|35|) as 

(f[H, h- *]) Mrs = ^tY. I dHdh M Rs[H, h, a\9]f[H, h; a] (36) 
M RS cr J 

where we defined the normalization constant M^ s = ^ a JdHdh Mng[H,h,tr\6]. 
Clearly ]im n ^oM^ s = 1. We use the above results to solve equation (|30|) for the 
functional distribution W[{P}], upon substituting the various RS expressions: 

A' 



M n RS 



f {dp} w[{p}] n p(a a , h a ) = yi I n { dj ? p ( j ?) wm}]} 

a = l fe>0 ' "* £=1 

x Y[ dH a d(<j a ,H a )e ih " [H «- e] Yl [Y, dhfP t (af,hf)e- iJ ^ h ^ +h "' T ^ 

a=\ J 1=1 af J 

n „ k 

{dp} n p(a a ,h a ) c -n^ c / n { dj i p ( j f) {UK wm}}) 

a = l k>0 ' J 1=1 

Z»[{J\,...,f\}] 

JdHd(a,H)e^ H ^ IlLiE^ fM t P t (<Tt, h e )e~ iJ ^+ h ^} 



l[s[p(a,h) 



Z[{Px,...,P k }] 

tH 
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where 

Z[{Pi, . . . , P k }} = I dHdh d ^ H)e ik ^ 

a J 

x J] {E / dh e P e (a e ,h e )e~ u ^+ h ^} (37) 

1=1 at J 

= 2tt YJ{ { E J&tPtfa, ht)e- iJ '~ h ">}d(a,J2 • 



a 1=1 " a e 

In the limit n — > both the normalization term Z n [{P±, . . . , Pk}] and the constant 
M^ s reduce to unity, and we find an equation for the functional distribution VF[{P}]: 

k 

' ' (38) 

t=i K 



W H P K = E ^ e " C / II {<VtP(Jt) {dPe} W[{P e }}} 
k>a ' J e=i 



]Js[p(a,h) 



Z[{P 1 ,...,P k }} 



In a similar fashion (see |Appcndix B| for details) we can compute the probability 
distributions D(s, h) and A[s, s'; h, h'; s] in RS ansatz. To compactify our formulae we 
define the Fourier transforms P{a\x) = Jdh P(a, h)e~ lhx , in terms of which we find 

D(s, h) = d(s, h) %z~ c J II {dJiPW {dPt}W[{P t }] } 



fc>0 



f=i 



Ut=i { Pt(<re\Jes)}6[h—J2e =1 Jecr e - 
X Z[{P 1) ...,P k }] 

A[s, h, ti- s] = \f~ C / II { A -hP(Ji) {dPt} W[{P £ }] } 



(39) 



k>0 



E C —f C U{ dJ r P ( J r){dQr}W[{Q r }]} 
m>0 ' r=l 

(jll^M^tlJe^nsw-^J^t-o-JsW,^ 

\e=i I a e ) 1=1 

m ( \ m 

111 ^2Qr(a- r \JrS) \ 5[h-^J' r (ir-e-Js'+2Js\d{s,h+2Js) 

=1 I oy J r=l 

k ( \ k 

EmE^I- 7 ^) \ 4 CT 'E Jtat + O+Ja') 



oo' £—1 \ (Tf> 



(40) 



Our theory requires solution of the saddle-point equations ()38ti39[) for the functional 
distribution VF[{P}] and the function d(s,h). These are functional relations, and 
is generally not possible to solve them analytically. Furthermore, the imaginary 
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arguments in (|38[) induce further complications in numerical solution. To simplify 
matters we assume that for x G 1R the Fourier transforms P(a\x) are real- valued, and 
we define a corresponding functional distribution 



W[{P}] = J {dP} W[{P}} Y[ 8 [P(o-\x) -fdh P(a, h)e- ih 



(41) 



W is normalized by construction, but the P{a\x) need not be. We transform our 
problem into the language of W by inserting (f3"5| into (|4"Tj) and integrating over {P}: 

W ^ = E P~ C I II {dJtP(Jt){dPz}W[{P e }]} (42) 



fc>0 



P(a\x) 



Z[{P 1 ,...,P k }} 



Our previous results (139M40P now take the form 



and 



D( S ,h) = d(s,h)J2^e- c [ Y[{dJ e P(J l ,){dP e }W[{P e }}} 

k>o ' J e=i 

nti {E CT , PtWts)} 6[h-j: k e=1 Ji°-i-0} 

Z[{P 1 ,...,P k }] 

k r & 

A[s,s';h,h';s\ = £ C -e~ c j J[ {dJ,P(J £ ) {dP e } W[{P e }}} 

L.>.n * ** n i 



(43) 



fc>0 



E h 6 ~ C / II { dJ rdP(J' r ) {dQr} W[{Qr}]} 
m>0 ' r=l 

(]l\^M^t\J^^)\s[h'-Y f J^t-0-JsW,h , ) 

\l=l [a t J 1=1 

m ( \ m 

Yl\ ^2Qr(a- r \JrS) \ S[h-^ f J' r a r -d-Js'+2Js\d(s,h+2Js) 

=1 I oy J r=l 

k ( \ k 

ELn E^i- 7 ^) \ d ( CT 'E Mi+o+ja') 



,(7(7* £ — 1 \ <j£ 



r—1 K. <r r 



(44) 



Equations (|42H44[) are the final analytic results within the replica symmetry theory. 
They complement and close the diffusion equation (|16p . We may now proceed to the 
solution of (|16p by iterating the following recipe from time t = onwards: at any 
time point t we use the instantaneous distribution D t (s, h) to solve equations (|42M3p 
numerically for W[{P}] and d(s, h) via a population dynamics algorithm 6], the result 
of which is then used to compute the kernel (|44|) . and to iterate (fl~6f over the next 
infinitesimal time step t — ► t + dt. 
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4. Application of the theory to the vertex cover problem 

In this section we show how the theory developed in the previous sections allows us to 
study analytically the Monte Carlo dynamics (extended with appropriate stochastic 
cooling schedules of the simulated annealing type) of the so-called minimum vertex 
cover optimization problem. 

4-1- The minimal vertex cover problem 

We start with the definitions. Let G = (V, E) be a graph defined by a set of iV vertices 
V = {1, . . . , N} and a set of undirected edges E = where i, j £ V and there is 

no distinction between and A vertex cover (VC) of a graph G is a subset 

Vvc Q V of vertices, such that for all edges (i,j) £ E either i £ Vvc or j £ Vvc 
or both. The vertices in Vvc are called covered, and those in V \ Vvc uncovered. 
Similarly, an edge (i,j) is said to be covered if at least one of the vertices {i,j} is in 
Vvc- The minimum vertex cover problem is the following optimization problem: find 
a vertex cover set Vvc of minimal cardinality, for a given graph G, and compute the 
fraction x c (G) = \Vvc\/N. The corresponding decision problem, to find whether for 
a given graph G a VC of fixed cardinality x = \Vyc\/N exists, is known to belong 
to the class of NP-complete problems [25]; i.e. it is conjectured that no algorithm of 
polynomial time complexity in N (or in the number of edges M) exists to solve it. All 
algorithms known to date indeed have exponential time complexity. The introduction 
of graph ensembles allows to study typical instances of the minimal VC problem, and 
quantify average properties. For instance, it was found that for large Poissonnian 
random graphs (J3J) the fraction of covered vertices in a minimal vertex cover x c (G) 
depends on the average connectivity c only, i.e. x c (G) = x c (c). Rigorous lower and 
upper bounds for x c (c) were derived in 26J: 

xi(c) < x c (c) < 1 - ln(c)/c (45) 

where xi(c) is a solution of 

xlnx+{l-x)la(l -x)- ^c(l - x) 2 = (46) 

The lower bound coincides with the annealed bound calculated within statistical 
mechanics, see e.g. [27] • The asymptotic form of x c (c) for large c was given in [2B] : 

x c (c) = 1 - 2[1 + ln(c) - ln(ln(c)) - ln(2)]/c + o{c~ l ) (47) 

Since in VC problems a vertex is either covered or uncovered, one can map the 
VC problem in to an Ising model log for any subset U C V we define a { = 1 
if i £ U, and <7j = — 1 if i £ U. We define a corresponding Hamiltonian for the state 
cr = (a%, . . . , tTjv) which simply counts the number of uncovered edges: 

iJ(o-) = ^c i ^ CT<) _ 1 ^.,_ 1 (48) 

i<j 

Solving the VC problem for a given relative cardinality x then reduces to minimizing 
H{cr) under the constraint JV = xN , which in terms of the Ising spins implies 

^I> = 2z-1 (49) 

i 

§ An alternative representation involves mapping the VC problem into a hard-sphere lattice gas 
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The study of VC has thereby been connected to the study of the ground states of 
an Ising spin system, within equilibrium statistical mechanics. This enabled the 
computation of the relative size of minimal VCs for typical large graph instances 
of FC ensembles, by averaging the free energy of the spin system over all graphs in 
the ensemble. Within the RS ansatz this resulted for Poisonnian graphs in [TQl [29] : 



«,(■:) (50, 

Here W(c) is the Lambert W- function [30], defined as the real solution of c = We w . 
This result implies that almost all graphs are coverable with xN vertices for x > x c (c) 
and not coverable for x < x c {c). It was proved to be exact for c < exp(l) [31] , but 
for c > exp(l) equation (|50[) underestimates the empirical values of x c (c) obtained 
by numeric simulation [10 , and for c > 20.7 it even violates the lower bound ([46]) . 
The explanation was found to be that at c = exp(l) the assumed RS breaks down 
PIJIHH], and replica symmetry breaking (RSB) occurs. The one-step RSB solution was 
computed in |32j via the cavity method; its agreement with the numeric results |10j 
improves upon the RS calculation (|50| and approaches ([47]) correctly for large c, yet 
the one-step RSB solution is still incorrect for c > exp(l) [32 . A more recent result 
obtained in [33] is in good agreement with both the numerical simulations [10| for 
c < 10.0, and with the asymptotic form (|47p for large c. The few and limited analytic 
studies of the dynamics of algorithmic solutions of the VC problem were carried out 
only for a simple backtracking algorithm |34j . and for more complex heuristic [35] 
algorithms. In this paper, in contrast, we consider a more physical dynamics, inspired 
by the connection with a ground state search in Ising spin systems. 



4-2. Dynamic replica analysis of the vertex cover problem 

In this section we analyze a Monte Carlo dynamics for the minimum VC problem 
of the type (O, where we allow the temperature T to vary over time such that 
limt^oo T(t) = 0, but sufficiently slowly so that in ([5]) we may simply substitute 
(3 — > (3(t). We map the VC problem into the Ising model (j48|) . and impose the 
constraint (|49|) in a 'soft' way, by adding an extra term to the Hamiltonian (f48|) 

H(cr) = yjcjj^, -i^,-! - A^(5 ffi! _i (51) 

i<j i 

(where A > 0) , which ensures that among the states cr that minimize (|5ip , those with 
the smallest sum X)i^o»,i (he. minimal cover) are preferred. The Glauber dynamics 
associated with the Hamiltonian ([ST]) is indeed of the type (pj, with the local field 

h i (<r) = J^c i j5 ajt _ 1 + d (52) 

where J = h and = — |A. The fields could also have been written in the more 
standard form hi{cr) = Jij&j + 0i, but this would have required site dependent 
random 9i which involve the connectivity variables c^. The consequences for our 
theory of changing the local fields from the conventional form @ to ([52"]) are minor. 
The diffusion equation (|16p remains unchanged; the only difference is in definition of 
the distribution (fTl"]) . which now becomes 

D{s,h;(r) = j^J2 Ss ^ h - J 12 Ci ^'- 1 -^ ( 53 ) 

i j^i 
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and the kernel (fTS)) , which changes to 

A[s,s';h,ti;s] = {^^S^S^Cij 5[h' -hi(a)]S[h-hj((T)-Js\^ (54) 

y 

where s £ {0,s'} and hi(er) is given by (|5^)) . Next we compute the consequences of 
defining (|52p within the replica calculations. This involves only minor alterations of 
the steps taken in section [31 and we readily obtain the new expression that replaces 
our previous l|28p (where in VC there is of course no longer a need to average over J): 

A[s,s';h,h';s]= (55) 
(^V/[^ffI''-5rJ^" iJE ' M ''.'- 1+l ' i '"'- ll ) Mji 

lim 



M,A/ 



The saddle-point equations ()29|30jl remain unaltered, with the (. . .)m averages given 
by equation (|25[) . but now the associated measure takes the new form 

M[H,h,cr\6] = exp Uh ■ [H - 9] - i^ D a {a a , H a ) (56) 

a 

+ cJ2 Uh P( ff / s fc , )[e- 1JE -^.^ +,l « , "«- 1 l-l]} 

The only changes to the earlier theory that are induced by the introduction of ([52]) 
are in the imaginary arguments of the exponential function in (|55|) and (|56p . We can 
therefore derive the RS version of the theory for VC dynamics simply by replacing 
o~a — * 3cr a ,-i an d P{Ji) —* 6(Je — J) in equations ()38M0p of section I3~2l This results in 

k 



w ti p yi - E \f~ c \ n {mywuPt}]} (57) 

fc>0 ' t=l 

n*[ i '( £r ' fc ) 



and, with the Fourier transforms = Jdh P(a, h)e~ lhx , 

D(s,h)= d( S ,h)J2^e- c [l[{{dP e }W[{Pe}}} (58) 

Ilti { £ g , gNJU}g - J£jU <^,-i - g] 

x Z[{p 1; ...,p fc }] 

fc>0 ' £=1 

E ^y e " c /n{w^[w^E^^i J ^-i)} ( 59 ) 

m>0 ' r=l cr r 

A; 

d[h'-J^a e ,-l-0-J6 s ^} d(s',ti) 



X 

m>0 
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m 

x 8[h— jy^ CTr) _i — 6 — J^ 3 / ; _i — Js] d(s, h - Js) 

r—l 

k k 

x J{{ y ^Pt{^\JSa,-i)]d{a,jY J ^ e ,-i+e+J5 al ^ 1 ) 



aa 1 £—1 <T£ 



77 i ill 



As before we may switch to a measure defined directly on the relevant Fourier 
transforms, which in the case of VC simplifies further due to the uniform bonds J: 



W[{P}}= J{dP}W[{P}]l[s[p(a\ 



JSa',-i) - / dh P(a,h)e-' 1 



hJ8„ 



(60) 



Our RS equations now acquire the following form: 



= E Y\ e ~° I II (61) 

fc>0 ' £=1 

_ r. ntl{E CT /K^I^,-l)}^^Etl^,-l+^+^',-l) 

x TT' 5 k^i^-'-i) L — = s 

D(s,/0 = d(«,fc)2^e-° /n{{dA}^[{^}]} 



(62) 



k>0 



1=1 



nti {£ CT4 Aw<Vi)}^ - ^Eti 



'cr/,-1 



and 



Z[{Pi,...,P k }] 

A[s,s';h,ti;s] = E^"fn{^M^E P '( ff «l*.-i)) ( 63 ) 

fc>0 ' ^ £=1 <xc 

x e ^T e_c /nlw^^^^E^^i^.- 1 )} 

m>0 r—l ov 

fe 

x J^£ CT , -i -6-J6 s -i] d(s',ti) 

1=1 

m 

x S[h - J)^^ f -i - - J5 s ' t -i - Js] d(s, h — Js) 

r=l 

fe fe 

E II { E A^iJ^-oJd^j^^-i+e+j^^) 



m m 

n{E < 3 r ^i j ^ , '- i )} d ( <r '' j E^-- i+e+j ^- 1 ) 



n -1 



Compared to the more general expression (|42|) , in the VC case (fBTjl the dimensionality 
of our problem has been reduced drastically as W is now a functional on the space of 
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2x2 matrices P(cr\ J8 a i Furthermore, the solutions of (|62l63p are of the following 
form, which is expected on physical grounds (given the non-random bonds J in VC): 

D(s,h) =Y^P{s,n)5{h-Jn-9) (64) 

n>0 

A[s,s';h,h';s]= ^ A[s, s'; n, n'} 6[ti -Jri -O-JS^t] 

n,n>0 

x 6[h-Jn-9-J6 s ,-i - Js] (65) 

where P(s,n) and A[s, s';n,n'} (with s,s' G { — 1,1} and n,n' <E {0,1,2,...}) are 
solved from 

p M= E^ c /n{{ d ^w^}]} ( 66 ) 

fc>o ' J e=i 

Ilti {£ g< PtWSs,-i)}d(s,Jn+e)s n ^ LiS<ri _ i 

and 

A[s,s';n,n'] = £ ^e~ c /j] {{m}W[{P e }\ £ PM (67) 

fc>0 ' £=1 cr* 

x e ^ e-c / n{i d ^^[W'-}]E^(^i j ^-i)} 

m>0 ' r=l ov 

x d( 5 ', Jn' + 0+J<J s ,_i) ^n',E? = i <V<.-i 
x d(s, Jn+6>+ J5 s ',_i) 8 ni j2^ =1 S ar ,-! 

k k 

x Ell [E^(^i J ^-i)] d ( CT ' J E^.-i+ 0+J ^'-i) 

_ (7 (7 1 £=1 <7t t=\ 

rn m - — 1 

X n [E ^Orl ^<r',-l)]d(ff' ^E S <rr,-l+0 + JS *,-l) 
r—1 cr r r—1 

The simplified form of the kernels (|64I65[) subsequently allows us to transform the 
main dynamical equation (Til))) , which is a PDE, into the following system of ordinary 
differential equations (see | Appendix C| for details): 

^P(s,0) = i [l+stanh[/30]]P(-s,O) - ^ [l-stanh[/?(9]] P(s, 0) 



+ -c [l+tsaih\pjri+p9+0J6 at _ 1 ]]A[s-l-,O,ri] 

n'>0 

- \-c Y [l-tsaih[f3M+0e+f3JS s -i}]A[s, 1; 0, n] (68) 



2 

n'>0 



whereas for n > we have 



— P(s,n) = i [l + stanh[/3Jn+/36»]]P(-s,n) - ~ [l-stanh[/3Jn+/30]] P(s, n) 
+ -c ^ [l-tanh[f3Jn'+f36+f3J6 s -i]]A[s,l;n-l,n'] 

n>>0 
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+ -c ^ [l + tanh[/3Jn'+/36»+/3J(5 s ,_i]]A[s-l ; n, n J 

n'>0 

- -c ^ [l-tanh[/3Jn'+/36i+/3J(5 s ,-i]]A[s,l;n,n'] 

?i'>0 

- 2 C 5Z [ 1 + tailll [/ 3J «'+/ 36, +/ 3J ^-i]]^[ s (69) 

n'>0 

Equations (|6"Tj) and (|66H69p are the final results of our dynamical replica analysis. They 
can be solved numerically, using population dynamics for the functional saddle-point 
equations (see |Appendix E| for details) and any standard method for the system of 
ordinary differential equations. If we allow for temperature adaptation, viz. (3 — > P(t), 
and restrict ourselves to those cooling protocols where (3(t) changes only on O(N ) 
time scales, we may simply make the replacement (3 — > f3(t) in the above equations. 

5. Tests of the VC theory against numerical simulations 

To test our theoretic predictions for the evolution of observables in the Glauber 
algorithm with stochastic cooling running on the VC problem, we compare the results 
of solving numerically the system of dynamical equations (|6T))) with the results of 
numerical simulations. We solve (|69p using a simple first-order Euler method, i.e. we 
iterate the iteration 

P i+l {s,n)=P l {s,n)+hT[s,n;P t {...);A t [...}], t t = Ih (70) 

where n € {0, 1, ... , L(c)} and T[. . .] is a short-hand for the right-hand side of (|68I69[) . 
Here L(c) denotes a suitable cut-off value that increases monotonically with the 
average connectivity c in (J5J), and < h -C 1. At each discrete time-step £ of this 
iteration we solve the RS equations (|61l66p via a population dynamics algorithm (see 
section [Appendix Ep and compute the kernel (|6T|) . Solving (|6T)|) requires a significant 
numerical effort, the bulk of which is devoted to solving equations (|61l66p . where 
the computation of a typical object (jE.ip requires typically 0(2 C ) basic operations. 
Although we expect that for sufficiently small h the changes in the statistical properties 
of the population between consecutive iterative time-steps in (|70p are small, the fact 
that the running time of the algorithm (|70p grows exponentially with c restricts 
the scope of simulation experiments. For each choice of control parameters our 
experimental protocol has been the following. First we generate a large random 
Poissonnian graph with the required connectivity c . Then we run the algorithm 
|T]) with the local fields (|52p , from an initial spin configuration where the individual 
spins are drawn randomly and independently from the distribution (|D. 1|) . We then let 
the system evolve according to the Glauber algorithm, but with the temperature T(t) 
decreasing in stages (to achieve stochastic cooling), while we record the evolution of 
two macroscopic order parameters, being the fraction x of covered vertices 

^^E-w ( 71 ) 

i 

and the energy density E : which is proportional to the fraction of uncovered edges, 

i<j 
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Figure 1. Left and middle: evolution of the fraction x and the energy density E in 
the VC algorithm with simulated annealing, for c = 0.5, J = 1.0 and 9 = —0.99. 
Time is measured in iterations per spin. Solid lines: RS theory. Dashed and 
dotted lines: average and average plus/minus standard deviation as measured 
over 100 simulation runs in systems with TV = 10 4 spins. The annealing schedule 
had four stages: (i) T = 2 for t G [0, 10], (ii) T = 1 for t 6 [10, 20], (iii) T = 0.5 for 
t 6 [20,30], (iv) T = for t e [30,40]. Right: histograms (RS theory) of the two 
field distributions P(±l,n) at f = 40, together with the corresponding simulation 
measurements (markers with error bars). 



For a state <r to represent an acceptable vertex cover it must have E(cr) = 0. For such 
a cover to be minimal we want in addition x(<r) to be as small as possible. 

In simulated annealing one starts a Monte Carlo dynamics at a high temperature 
T(0), and then lowers it slowly in stages, allowing the system to equilibrate effectively 
along the way. The objective is for the algorithm not to get stuck in states that are 
only locally but not globally optimal. Determining the best cooling schedule T(t) 
for achieving this, however, is highly nontrivial; furthermore, equilibration times in 
VC-type optimization problems can scale exponentially in the system size. Here we 
did not attempt to optimize the cooling protocol but focused on the VC dynamics for 
simple step-wise temperature reductions. Our numerical simulations where carried 
out on random graphs with N = 10, 000 vertices, with average connectivities c £ 
{0.5,1,1.5,2,2.5,3,3.5}. The values of the parameters J = 1, 9 = —0.99, the initial 
covered fraction x = N^ 1 <5 CTi (o),i = 0.9 and the temperatures T G {2,1,0.5,0} 
(reduced in steps) were identical in all simulations. The initial conditions for (f69|) and 
the population dynamics were computed via equations (|D.3HD.5j) of |Appcndix D| The 
size of the population was M = 10, 000 and the number of iterations typically needed 
for the population dynamics to converge was of order 107V. 

In figures [TJH we compare the data obtained in our numerical simulations for 
c e {0.5,1,2,3} with the results of solving (|69|) numerically. We observe that the 
overall agreement between theory and simulations is excellent. The RS theory also 
predicts correctly the joint spin-field statistics P t (s,n), see the right panels in figures 
QUI note the different vertical scales. The deviations between theory and simulations 
are (as usual in DRT) confined to intermediate times, and limited to low temperatures 
in combination with high average connectivity, as shown in the insets of figures [3] and 
[4] but even there they remain within the error bars of the simulation data. Finally, in 
figure [5j we compare our results for the fraction of covered vertices x as measured at 
termination of the algorithm with the result (|50|) of equilibrium statistical mechanics, 
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P(l,n) 



Figure 2. Left and middle: evolution of the fraction x and the energy density E 
in the VC algorithm with simulated annealing, for c = 1, J = 1.0 and 9 = —0.99. 
Time is measured in iterations per spin. Solid lines: RS theory. Dashed and 
dotted lines: average and average plus/minus standard deviation as measured 
over 100 simulation runs in systems with TV = 10 4 spins. The annealing schedule 
had four stages: (i) T = 2 for t g [0, 10], (ii) T = 1 for t g [10, 20], (iii) T = 0.5 for 
t g [20,40], (iv) T = for t g [40,50]. Right: histograms (RS theory) of the two 
field distributions P(±l,n) at f = 50, together with the corresponding simulation 
measurements (markers with error bars). 




Figure 3. Left and middle: evolution of the fraction x and the energy density E 
in the VC algorithm with simulated annealing, for c = 2, J = 1.0 and 9 = —0.99. 
Time is measured in iterations per spin. Solid lines: RS theory. Dashed and 
dotted lines: average and average plus/minus standard deviation as measured 
over 100 simulation runs in systems with TV = 10 4 spins. The annealing schedule 
had four stages: (i) T = 2 for t g [0, 10], (ii) T = 1 for t g [10,20], (iii) T = 0.5 
for t g [20, 50], (iv) T = for t g [50,60]. The inset in the left figure shows an 
enlargement of the region t g [20,30], where the largest deviation between theory 
and simulation for x is observed. Right: histograms (RS theory) of the two field 
distributions P(±l,n) at t = 60, together with the corresponding simulation 
measurements (markers with error bars). 



as obtained (within the replica-symmetry ansatz) in 10J . Our data (markers, with 
virtually no difference between the observed values in simulations and the prediction 
of our dynamical theory) are seen to be close to the equilibrium prediction (|50| (solid 
curve), but they overestimate slightly the size of minimal vertex covers. This type 
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P(-l,n) 
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Figure 4. Left and middle: evolution of the fraction x and the energy density E 
in the VC algorithm with simulated annealing, for c = 3, J = 1.0 and 9 = —0.99. 
Time is measured in iterations per spin. Solid lines: RS theory. Dashed and 
dotted lines: average and average plus/minus standard deviation as measured 
over 100 simulation runs in systems with TV = 10 4 spins. The annealing schedule 
had four stages: (i) T = 2 for t £ [0, 10], (ii) T = 1 for t g [10,20], (iii) T = 0.5 
for t e [20, 50], (iv) T = for t G [50,60]. The inset in the left figure shows an 
enlargement of the region t £ [20,30], where the largest deviation between theory 
and simulation for x is observed. Right: histograms (RS theory) of the two field 
distributions P(±l,n) at t = 60, together with the corresponding simulation 
measurements (markers with error bars). 




Figure 5. The fraction x(c) of covered vertices in a minimal vertex cover as a 
function of the average connectivity c. Solid line: prediction of a static replica- 
symmetric calculation (exact for c < e; the value c = e is shown as a vertical 
dashed line), as obtained in |10| . Symbols: the predicted final fraction x of 
covered vertices in the vertex cover obtained from the present dynamics (when 
we have arrived at T = 0), according to the RS dynamical replica method (which 
agrees perfectly with the simulations). 



of behaviour is not unusual in simulations, and suggests that more sophisticated 
annealing schemes must be used to achieve equilibration. Slightly more unexpected is 
the fact that our RS theory exhibits a similar overestimation of x. This, in combination 
with the fact that the static RS equations of 10J can be shown to constitute a 
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stationary solution of our present dynamical replica equations, suggests (at least to 
the left of the vertical dashed line, where replica symmetry should hold) that the 
time required for true equilibration diverges with N. Both the static RS replica 
equations and the long time limit of the dynamical RS replica equations represent 
distinct stationary solutions of the dynamical formalism, and the observed differences 
in x{c) are manifestations of the non-commuting of the limits N — > oo and t — > oo. 

6. Discussion 

In this paper we have studied the sequential dynamics of finitely connected Ising 
spin models with random bonds, on Poissonnian random graphs. Starting from the 
microscopic master equation we derived a dynamic equation for the joint spin-field 
probability distribution, which is exact in the infinite system size limit, but not closed. 
We then followed the usual prescriptions and assumptions of dynamic replica theory 
[2"2"] in order to close this equation. The result is a set of nontrivial coupled diffusion 
equations, in which the evaluation of the driving forces requires the solution of a 
saddle-point problem at each instance of time. The latter saddle-point equations are 
of a functional nature, and are derived within the replica-symmetric (RS) ansatz; they 
can be solved numerically by a conventional population dynamics algorithm [6] . 

As a first application, we have applied our dynamical theory to the dynamics of 
a simulated annealing algorithm (Glauber-type dynamics with a stepwise stochastic 
cooling schedule) when running to find a solution of the so-called minimal vertex 
cover (VC) problem [TO], on finitely connected Poissonnian random graphs. In this 
problem the local fields are essentially integer- valued, which simplifies our dynamical 
equations. We have derived dynamic equation for the joint probability of spins and 
nonnegative integer fields. Upon solving the equations of our theory numerically and 
comparing the results with the outcome of numerical simulations of the algorithm, we 
find excellent agreement between theory and experiment. 

When compared to e.g. the generating functional analysis method (GFA), the 
advantage of dynamical replica theory (DRT) is that, unlike GFA, it does not give 
an effective number of scalar order parameters that grows exponentially with time. 
Although also in its present forrrjj], with the joint spin-field distribution as the core 
dynamical order parameter, the DRT method is not exact, it is certainly much more 
accurate than e.g. any simple two-parameter theory [17] . and it can be systematically 
improved further by increasing the order parameter set [T5], although at a numerical 
cost. We believe its wide applicability to be the main advantage of the dynamical 
theory presented in this paper. The formalism can be extended relatively easily to 
include, for instance, directed or non-Poissonnian random graphs. 
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|| For parallel stochastic dynamics it can be shown that an exact formulation of DRT is possible 
for any discrete spin model, with an effective number of scalar order parameters that grows at most 
linearly with time 36 . 
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Appendix A. Averaging over disorder 

In this section we give the details of the disorder averaging in equation (|20[) that brings 
us to equation (f2~Tj) . Firstly, we rewrite slightly the term within the angular brackets, 
exploiting the symmetry of CijJij under index permutations i <-> j: 

(■ ■ ■) {CtJ . hj} = (ci2S[h - Hi + 2J 12 s}e-^ h ? h ^) {CijJtj} 

= e-wS-«^ {ci2 S[h-H 1 2 +2J 12 S}e- i ^^ J ^^ h f^) {CijJij} (A.l) 

= e - w ^ h ° (c 12 6[h-H 1 2 +2J 12 S]e- i ^<> <* J » £J*?"?+*?<'?]) {cy Jy} 

We then average over the connectivity disorder {cj, ■}, which is defined by followed 
by the bond disorder { Jij}: 

(...) {cj . /i . } = ^e- ie S«^?(^[/i-ffl+2J 12 s]e- Ul2 ^^ CT =+^>"] 
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= -^e- w ^«ih? J dJ p(j) S[h-H^ + 2Js}e- u ^ [h °^ + '^ a ° ] 
x II {|/dJP(J)e-^^Kl + l-|} (A.2) 

Finally, we re-exponentiate the last line of the above expression, giving 

(. . .} {c . . Jij} = -le""'^'* / dJ P(J) ( 5[/i-i7 2 1 + 2JS]e- iJ ^=[' ! ° <T ? + ' l ? <T fl 



x exp 



JLY,lJ d JP ( J)e-^>! 



+ 0(1) (A.3) 



Appendix B. Calculation of the RS saddle-point equations 



We compute the RS versions of the kernel ([28)) and the saddle-point equation 
Assuming replica-symmetry transforms in these equations the averages over the 
effective measure, (. . .)m — * (■ ■ ^Mrsj with the definition (|55|) . In (|28I29[) this gives 

D(s, h) = j±rJ2 ^ e ' c [ II { d Jt p (Jt.) m} w[{P e }}} 



k>a 



i=i 



dffid/ii^cri, H^^-^Ss^Sih - Hi] 

k 

n {e / d^d^dt^,^)^''-^- 9 ' 

x 11 [E / d W(<^ ^)e- iJ '[ & -^+^] } 



1 cr 
k 



k>0 



E |t e " c / II { dJ z p W {m} w[{p e }} } 



Z[{p 1; ...,p fe }]' 



M 



IIS 



k „ 

Ell [E MiPfaMe 



a 1=1 
k 



d(a, Ji(Ti+0)S Sia 8[h-Y^ 



k>0 



E |f e_C / II {*Jt p {Jt) {dPi} W[{P,}} } 



Z[{P u ...,P k }f- 1 



M 



RS 



rf(s»n[E [ d h e P e (a e ,h e )e- iJ ^S[h~J2j^i-d} (B.l) 



-1 o~t 
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and 



A[s, s'; h, h ; s] 



A[s, s'; h, h'; s] 



(B.2) 



where 



X E^ e ~ 

k>0 



J2 aa , JdHdH' A[a,a';H,H';s] 

A[ 8 ,J;h,ti;a\ = (6 s ,, ai 6 s ^S[h'-H 1 ]6{h-H{+2Js}e- iJ ^ <T ' +h ' ^) j, Mrs ,m' rs 

= TtW / dJP ( J ) E [dHdH'dhdh8 sl ^ 1 5 s ^ i 8[ti-H 1 }5[h-H' 1 + 2Js] 
M RS J a-cr, J 

k 

Y[{dJiP(Jt){dP t }W[{Pi}]} 

e=i 

a k „ 

Yl d(<T a ,H a )e ih ^ H «- 9 - Ja '^ J] [E / dhfP e [af,hf)e- iJ ^ haa " + " h "' T ^ 

a=l 1=1 cr° 

E S°~ C l A { dJ rP(Jr) {dQr} W[{Q r }}} 
m>0 m ' r=l 

n m p 

[] d(*' a ,H' a )j^ H '<*- 9 - J ^ I] [E / dh?Q r (a?rK)e- iJ ^>? +h >'<> 
a=l r—1 of J 

^ E ^ e ~ c / ft {v t p(Ji) {m} w[{ Pi }}} 

RS k>0 ' J 1=1 

E —\ e ~ C ff[{dJrP(Jr){dQ r }W[{Q r }}} I dJ P(J) 
rn>0 m ' J r=l •* 



d(a,^2j e a e + 8 + Ja') 



Ell [E <&tPt{<Tt,'hi)e 

rn p m 

II [E dh r Q r {a r ,h r )e- Urhra ' d(a' , ^ J r a r +9+ J a) 



r—1 <r r 



5 s , a >5[h-^J r <j r -6-Jo-+2Jd\8 sl ^8[ti-^J l <j t -6-J<j'\ (B.3) 



i=i 



In the replica limit mOwe get lining M^ s = lim^o Z[{Pi, . . . , Pk}] n — 1, and 
also lim„^o Ylaa' JdHdH' A[a, a' ; H, H'; s] = 1. As a result, the above expressions 
(|B.1IB.2|) reduce to equations (|39|40p . 



Appendix C. Reduction of PDE to the system of ODEs 

In this section we show how the diffusion equation (|16|) , written in terms of the kernels 
(|53l54p . can be reduced to a system of ordinary differential equations. The discrete 
nature of the fields (|52|) allows to write (|53I54[) in terms of the probability distributions 
(|64jl and (|65| . Inserting (|64|) and (|65|) into both sides of (fl6|) gives 

d 1 

— 2_,P t [s,n)8(h-Jn-6) = - [l + stanh[(3h}}^2 P t (-s,n)8(h- Jn-0) 

n>0 n>0 
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- - [l-st&nh[f3h]]'^2Pt(s,n)S(h-Jn-e) 



n>Q 



+ \ C Y1 J dh ' [l-s'tanh[/3/i']] Ms, 

s ; n, n J 

s' nn'>0 

x 6[h'-J(n'+5 Si _ 1 )-d]S[h-J(n+S sli _ 1 +s')-9} 
J dh ' [l-s'twh[f3ti]] A[s,s';n,n'} 

s' nn' >0 

x S[ti-J(n'+S s ,-i)-d]8[h-J(n+S s > (C.l) 

In the left-hand side we move the time derivative inside the summation, and in the 
right-hand side we sum over s' and integrate over h! . This leads to 

^Pt(s,n)5(h-Jn-6) = S(h- Jn-6)\^ [1+s tanh[ ( 5( Jn+0)]] P t {-s,n) 

i [l-stanh[/3(Jn+0)]]P t (s,n)j 
1 -cY / S(h-Jn-9)J2{ 



n>0 



2 

n>0 n'>0 



[l+tanh[/3(Jn / +0+ J5 S _i)]]A t [s,-l; n, n'] 
- [l-tanh[/3(Jn' + 6» + J(5 S! _ 1 )]]A t [s,l;n,n']| 

n>l n'>0 

[l-tanh[/3(Jn' + 6»+J^ S: _i)]]A t [s,l;n-l,?i'] 
- [l+tanh[/3(Jn'+e+J^ s _x)]]At[s,-l;n-l,n']} (C.2) 
It follows from the above that the evolution in time of Pt(s,n) is governed by (f69|) . 



Appendix D. Initial conditions 

In this section we compute the values of the probability distribution P(s,n), the 
functional distribution W^[{P}] and the function d(s, Jn + 9) at time t = 0. We choose 
an initial state of the system in which all individual spin values are drawn randomly 
and independently, according to 

N 1 1 
Po(tr) = J[{ 2 (1 + m °^ 1 + 2 (1 " m °^-^ ( D - X ) 

i=l 

where mo G [— 1,1] is the prescribed initial magnetization. It follows that the joint 
spin-field distribution (f53|) at t = is given by 

1 N 

rr i 
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-(1 + smo) 2^ -j exp 



(1-mo) <S[/i-Jn-0] (D.2) 



n>0 



The initial conditions for the system (|69|) follow directly from the above expression: 

1.. J4c(l-m )]" , 1 



P (s,«) = ^(1 + sm )-^ : — — exp[-~c(l - m )] 



(D.3) 



2 V " n\ 2 

Furthermore, we see that the joint spin-field distribution (|D.2p indeed takes the desired 
form of the saddle-point equation (|62p , with the functional distribution 



W[{P}] = [] S[P(a\ J5 &; -i) - ~(1 + ™ ) 



and with 



<i(s, Jn + f 



1 



-(1 + sm ) 



(DA) 



(D.5) 



It is a trivial matter to show that (|D.4ID.5[) are indeed the solutions of equation ([61 



Appendix E. Population dynamics 



The functional saddle-point equations (|61|66|) cannot in general general be solved 
analytically (one trivial exception is the infinite temperature regime). We therefore 
resort to the so-called population dynamics algorithm [6] to obtain solutions 
numerically, solving equations (|61[) and (|66[) simultaneously for the functional 
distribution W"[{P}] and the function d(s, Jn + 9), given the (known) values of the 
probability distribution Pt(s, n) at time t. We create a population of A/" 2x2 matrices 
Pi(a\J8 a i ,-\), where i = 1 . . .A/", and we initialize the numbers d(s, Jn + 9), where 
s E { — 1,1} and ri 6 {0,1,2,...}. We then execute an iterative process, whereby 
at each step we update the population of matrices and the numbers d(s, Jn + 9) as 
follows: 

(i) a number k is drawn from the Poisson distribution P c (k) ([3]) 

(ii) k members P{(a\ JS a i ^i) are selected randomly and independently from the 
population 

(hi) a new value for P(a\JS a ' is calculated according to 



,{a\J8 a 















r e\Jo~a-i) 





(E.l) 



(iv) a member of the population is selected randomly, and replaced with the newly 
computed value P new (a\JSa> 

(v) a new function d(s, Jn + 9) is computed according to 

E c ^ c I ri {{m}w[{Pt}}} (e.2) 



,(s, Jn + 9) = P t (s, n) x 



fe>0 





= 1 [Eff; Pi(<Ti\J8 a ,- 


-0 






"E„A(^|J^,-i)" 





Analysis of processes on finitely connected graphs I: vertex covering 



20 



Here averaging over the functional measure W is denned as averaging over the actual 
instantaneous population of 2x2 matrices. This iteration is repeated until the values of 
the function d(s, Jn + 0) and the statistical properties of the population are stationary. 
The population measure W will now be an estimate of the functional distribution (I61|) , 
and the function d(s, Jn + 6) is a hxed point of the iteration equation (|E.2[) . i.e. a 
solution of our original saddle-point equation. 



